1/4 de la note finale est liée à la mise en forme :
Ce TP reprend l'exemple d'un médecin et de ses vaccins. Vous allez comparer plusieurs stratégies et trouver celle optimale. Un TP se fait seul ou en binôme. Aucun groupe de plus de 2 personnes.
Vous allez rendre le TP depuis un lien GitHub avec ce notebook mais une version du rapport exportée en PDF & HTML.
! pip install matplotlib tqdm numpy ipympl opencv-python
!jupyter labextension install @jupyter-widgets/jupyterlab-manager
!jupyter labextension install jupyter-matplotlib%load_ext autoreload
%autoreload 2
%matplotlib inline
import typing as t
import math
import torch
import numpy as np
from tqdm.auto import trange, tqdm
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.animation as animation
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
import cv2
from IPython.display import display, clear_output
torch.random.manual_seed(0)
K = 5 # num armsclass ArmBernoulli:
def __init__(self, p: float):
"""
Vaccine treatment following a Bernoulli law (mean is p and variance is p(1-p)
Args:
p (float): mean parameter
>>> torch.random.manual_seed(random_state)
>>> arm = ArmBernoulli(0.5)
>>> arm.sample(5)
tensor([ True, False, True, True, True])
"""
self.immunity_rate = p
def sample(self, n: int = 1):
return torch.rand(n) < self.immunity_rate
def __repr__(self):
return f'<ArmBernoulli p={self.immunity_rate}'
def generate_arms(num_arms: int):
means = torch.rand(K)
MAB = [ArmBernoulli(m) for m in means]
assert MAB[0].immunity_rate == means[0]
assert (MAB[0].sample(10) <= 1).all() and (MAB[0].sample(10) >= 0).all()
return MAB
MAB = generate_arms(K)Ce TP reprend l'exemple du médecin présenté en cours.
Q1. Créez une fonction pour trouver μ* à partir d'un
MAB. Comment est définie la récompense Rk ? Que
représente concrètement le regret dans le contexte de ce TP
?
#set random_state to 0
torch.manual_seed(0)
MAB = generate_arms(K)
best_arm = max(MAB, key=lambda MAB: MAB.immunity_rate)
arms_immunity_rate = [arm.immunity_rate.item() for arm in MAB]
print(f'Immunty_rate for random_state == 0 is : {arms_immunity_rate}')
print('----')
print(f'Best vaccin is : {best_arm}')Immunty_rate for random_state == 0 is : [0.49625658988952637, 0.7682217955589294, 0.08847743272781372, 0.13203048706054688, 0.30742281675338745]
----
Best vaccin is : <ArmBernoulli p=0.7682217955589294
Définition récompense R_k
La récompense R_k est définit par une loi de Bernouilli.
Représentation du regret
Le regret représente la différence de personne qui aurait du être
immunisé après l'injection du vaccin mais qui aurait du l'être
Note importante : pour la suite, les résultats
seront généralement réalisés avec 100 initialisations différentes du MAB
(tous les MAB ont 5 vaccins mais des taux d'immunistation différent)
pour réduire le bruit de simulation. Concrètement, on exécutera au moins
100x generate_arms.
Le médecin fonctionne sur deux phases :
$$\hat{\mu_i}[0\rightarrow N] = \frac{1}{T_i} \sum_{k=0}^{N-1} \chi_{v_k,i}R_k,$$
avec $T_i = \sum_{k=0}^{N-1} \chi_{v_k,i}$.
Q2. Implémentez cette solution avec N = 50 et M = 500 et testez-la avec 100 MAB. On souhaite savoir si vous trouvez le vaccin optimal à l'issue d'une phase d'exploration. Quelle est l'espérance empirique de cette variable ? Et son écart-type ? Calculez de même l'espérance et l'écart-type du regret sur vos 100 simulations.
Pour rappel, le regret est défini par :
$$r_n = n\mu^* - \sum_{k=0}^{n-1} R_k$$
Attention : n est le nombre total de patients, donc ici N + M.
#Déclaration des constantes
N = 50
M = 500
number_MAB = 100def getBestVaccin(MAB):
best_arm = max(MAB, key=lambda MAB: MAB.immunity_rate)
best_index = MAB.index(best_arm)
return best_index
def getBestResult(vac_results):
return vac_results.index(max(vac_results))
def Vaccin_arms(MAB, n):
patient = []
for i in range(len(MAB)):
vaccins = MAB[i].sample(n)
vaccins = [1 if vaccin else 0 for vaccin in vaccins]
patient += (vaccins)
return patient
def Exploration(MAB, n=N): #return [1,0,...],
Ti = n // len(MAB) #Distribution uniforme
vaccins = Vaccin_arms(MAB, Ti)
vac_results = []
for i in range(len(MAB)):
sum = 0
for j in range(Ti):
ind = i * Ti + j
sum += vaccins[ind]
vac_results.append(sum / Ti)
return vac_results
def Exploitation(MAB, vac_results, m=M): #return [1,0 ...], for best vaccin
arg_max = getBestResult(vac_results)
vaccins = MAB[arg_max].sample(m)
vaccins = [1 if vaccin else 0 for vaccin in vaccins]
return vaccins
def gloutonne(number_MAB=100):
exploration_var = []
regret = []
for _ in range(number_MAB):
MAB = generate_arms(K)
exploration = Exploration(MAB, N)
best_arm = getBestVaccin(MAB)
if (best_arm == getBestResult(exploration)):
exploration_var.append(1)
else:
exploration_var.append(0)
exploitation = Exploitation(MAB, exploration, M)
Ti = N // len(MAB)
r_N = (N + M) * MAB[best_arm].immunity_rate - ((exploration[best_arm] * Ti) + sum(exploitation))
regret.append(r_N)
return exploration_var, regret
exploration_var, regret = gloutonne(number_MAB)
print(f'Espérence empirique de l\'exploration : {np.mean(exploration_var)}')
print(f'Ecart-type de l\'exploration : {np.std(exploration_var)}')
print(f'Espérence empirique du regret : {np.mean(regret)}')
print(f'Ecart-type du regret : {np.std(regret)}')Espérence empirique de l'exploration : 0.75
Ecart-type de l'exploration : 0.4330127018922193
Espérence empirique du regret : 45.385528564453125
Ecart-type du regret : 34.105445861816406
Q3. On étudie maintenant l'influence de la taille du training set N. On considère que N+M est une constante, puis on fait varier N entre K et M. Calculez le regret pour ces différentes tailles du training set différents MAB et representez le regret moyen, le regret min et max (vous devriez trouver une courbe en U ou en V pour le regret moyen). Quelle est la taille optimale du training set ?
# Constantes
n_total = N + M
def gloutonne_2():
regret = []
n = K #5
m = n_total - n #545
MAB = generate_arms(K)
while n <= m:
exploration = Exploration(MAB, n)
exploitation = Exploitation(MAB, exploration, m)
#regret
Ti = n // len(MAB)
best_arm = getBestVaccin(MAB)
r_N = (n + m) * MAB[best_arm].immunity_rate - ((exploration[best_arm] * Ti) + sum(exploitation))
regret.append(r_N.item())
n = n + K
m = m - K
return regret
import matplotlib.pyplot as plt
def plot_regret(regret):
means = [np.mean(regret)] * len(regret)
plt.figure(figsize=(10, 6))
plt.plot(np.arange(0, len(regret)*K, K), regret, marker='o', linestyle='-', color='b')
plt.plot(np.arange(0, len(regret)*K, K), means, linestyle='-', color='r')
plt.xlabel('Valeur de N')
plt.ylabel('Regret')
plt.title('Courbe du Regret')
plt.grid(True)
plt.show()regret_v2 = gloutonne_2()plot_regret(regret_v2)
print('regret min :', min(regret_v2), ' | regret max :', max(regret_v2), ' | regret moyen :', np.mean(regret_v2))
print('La taille optimal du training set N est :', np.argmin(regret_v2) * K)
regret min : -5.828857421875 | regret max : 206.171142578125 | regret moyen : 101.51659712357954
La taille optimal du training set N est : 0
Q4. On propose d'améliorer l'algorithme précédant en mettant à jour les taux d'immunisation empiriques μ̂i pendant la phase d'exploitation (algorithme greedy). Concrètement, à chaque nouvel patient, on lui administre le meilleur vaccin selon les stats. Notez vous une amélioration du regret ? Proposez un exemple où les taux d'immunisation du MAB ne changent rien.
def Exploitation_greedy(MAB, immunity_rate):
vaccins = []
Ti = N // len(MAB)
hist = [Ti] * len(MAB) # size -> 5
for _ in range(M):
best_arm = getBestResult(immunity_rate)
arm = MAB[best_arm].sample(1).item()
hist[best_arm] += 1
immunity_rate[best_arm] = (immunity_rate[best_arm] * (hist[best_arm] - 1) + arm) / hist[best_arm]
vaccins.append(arm)
vaccins = [1 if vaccin else 0 for vaccin in vaccins]
return vaccins
def gloutonne_greedy(number_MAB=100):
regret = []
for _ in range(number_MAB):
MAB = generate_arms(5)
exploration = Exploration(MAB)
best_arm = getBestVaccin(MAB)
exploitation = Exploitation_greedy(MAB, exploration)
Ti = N // len(MAB)
r_N = (N + M) * MAB[best_arm].immunity_rate - ((exploration[best_arm] * Ti) + sum(exploitation))
regret.append(r_N)
return regret
def plot_regret_greedy(regret):
means = [np.mean(regret)] * len(regret)
plt.figure(figsize=(10, 6))
plt.plot(np.arange(0, len(regret), 1), regret, marker='o', linestyle='-', color='b')
plt.plot(np.arange(0, len(regret),1), means, linestyle='-', color='r')
plt.xlabel('Itération')
plt.ylabel('Regret Moyen')
plt.title('Courbe du Regret')
plt.grid(True)
plt.show()
regret_greedy = gloutonne_greedy(number_MAB)
plot_regret_greedy(regret_greedy)
print(f'Espérence empirique du regret : {np.mean(regret_greedy)}')
print(f'Ecart-type du regret : {np.std(regret_greedy)}')
Espérence empirique du regret : 38.36371994018555
Ecart-type du regret : 13.835759162902832
Réponse
La figure montre le résultats sur 100 itérations contenant tous 5 MAB
aléatoires.
Nous remarquons que l'écart-type du regret diminue drastiquement.
Quant à la moyenne du regret, elle ne diminue que légèrement (< 10).
Les résultats sont donc dans l'ensemble meilleur avec cette algorithme
(greedy).
Cette mise à jour ne change rien dans le cadre où le meilleur vaccin théorique est trouvé dès le début (à la première itération).
Q5. Nouvelle amélioration : à chaque nouveau patient, on choisit si on lui administre le meilleur vaccin avec une probabilité ϵ ou un vaccin aléatoire (p = 1 − ϵ). Vérifiez si vous obtenez un meilleur résultat avec N = 0 ou N > 0. À votre avis, à quoi sert ϵ ?
#Constantes
eps = 0.9
def Exploitation_proba(MAB, immunity_rate, n=N, m=M):
vaccins = []
Ti = n // len(MAB)
hist = [Ti] * len(MAB)
for _ in range(m):
if (torch.rand(1) > eps):
best_arm = torch.randint(0, len(MAB), (1,)).item()
else:
best_arm = getBestResult(immunity_rate)
arm = MAB[best_arm].sample(1).item()
hist[best_arm] += 1
immunity_rate[best_arm] = (immunity_rate[best_arm] * (hist[best_arm] - 1) + arm) / hist[best_arm]
vaccins.append(arm)
vaccins = [1 if vaccin else 0 for vaccin in vaccins]
return vaccins
def gloutonne_proba(number_MAB=100, n=N, m=M):
regret = []
for _ in range(number_MAB):
MAB = generate_arms(K)
if (n == 0):
exploration = [0] * len(MAB)
else:
exploration = Exploration(MAB, n)
best_arm = getBestVaccin(MAB)
arg_max = getBestResult(exploration)
exploitation = Exploitation_proba(MAB, exploration, n=n)
Ti = n // len(MAB)
r_N = (n + m) * MAB[best_arm].immunity_rate - ((exploration[arg_max] * Ti) + sum(exploitation))
regret.append(r_N)
return regret
regret_eps= gloutonne_proba(number_MAB, n=N)
regret_eps_0 = gloutonne_proba(number_MAB, n=0)
print(f'Moyenne du regret pour N > 0 : {np.mean(regret_eps)}')
print(f'Moyenne du regret pour N = 0 : {np.mean(regret_eps_0)}')Moyenne du regret pour N > 0 : 54.599342346191406
Moyenne du regret pour N = 0 : 36.695518493652344
Réponse
Nous obtenons de meilleurs résulats avec N = 0. Cependant, le regret est moins bon avec l'amélioration qu'avec la méthode greedy classique
Epsilon permet de tester des d'autres vaccins, afin de ne pas être bloqué sur un maximum local (recherche ailleurs), c'est à dire, il établie une recherche aléatoire pour trouver le le bras avec la meilleure immunity rate.
Lai et Robbins [Lai et Robbins, 1985] considère une classe d'algorithmes π pour résoudre ce type de problèmes.
Ils ont trouvé une borne inférieure sur les récompenses cumulées en valeur asymptotique :
$$\lim_{n\rightarrow \infty} \inf_{\pi} \frac{\sum_{k=0}^{n-1} R_k}{\log n} \geq \sum_{i~\text{tel que}~\mu_i \lt \mu^*} \frac{\mu^∗−\mu_i}{\text{KL}(\mu_i, \mu^*)} :=C(\mu)$$
avec KL(x,y) = xlog (x/y) + (1−x)log ((1−x)/(1−y)) (distance de Kullback-Leibler) et $\sum_{k=0}^{n-1} R_k$ la récompense obtenue sur n patients.
Q6. Justifiez pourquoi on peut en déduire que le regret d'un algorithme raisonnable sera au pire logarithmique.
Réponse
La borne inférieure de Lai et Robbins montre que la récompense cumulée
d'un algorithme efficace suit asymptotiquement une croissance
proportionnelle à log(n)
:
$$ \lim_{n \to \infty} \inf_{\pi} \frac{\sum_{k=0}^{n-1} R_k}{\log n} \geq C(\mu) $$
Le regret R(n) est défini comme la différence entre la récompense maximale que l'on pourrait obtenir en jouant toujours le meilleur bras (celui avec la plus grande espérance, μ*) et la récompense effectivement obtenue par l'algorithme. Autrement dit :
$$R(n) = n\mu^* - \mathbb{E}\left[\sum_{k=0}^{n-1} R_k\right]$$
En utilisant la borne inférieure sur la somme des récompenses obtenues $\sum_{k=0}^{n-1} R_k$, on peut en déduire que :
R(n) ≥ nμ* − C(μ)log n
Cela signifie que, pour un algorithme raisonnable, le regret ne peut pas croître plus vite qu'une fonction logarithmique en n, c'est-à-dire que R(n) = O(logn) dans le pire des cas. En d'autres termes, tout algorithme efficace dans ce cadre aura un regret au plus logarithmique.
Ainsi, la borne inférieure de Lai et Robbins, en montrant que la récompense cumulée doit croître au moins comme C(μ)log n, implique que le regret d'un algorithme raisonnable sera également logarithmique dans n. Un algorithme "raisonnable" ne pourra donc pas avoir un regret asymptotiquement pire que O(logn).
Q7. Tracez le regret issu de la borne de Lai & Robbins et comparez le au regret obtenu avec l'algorithme glouton.
#Lay & Robbins
def kl_divergence(x, y):
return x * np.log(x / y) + (1 - x) * np.log((1 - x) / (1 - y))
def lower_bound(n, immunity_rate, mu_star):
C_mu = 0
for mu_i in immunity_rate:
if mu_i < mu_star:
C_mu += (mu_star - mu_i) / kl_divergence(mu_i, mu_star)
return C_mu * np.log(n)
def lai_robbins_regret(n, C_mu, mu_star):
return n * mu_star - C_mu
def lai_robbins(MAB, n):
immunity_rate = [arm.immunity_rate.item() for arm in MAB]
mu_star = max(immunity_rate)
C_mu = lower_bound(n, immunity_rate, mu_star)
regret = lai_robbins_regret(n+1, C_mu, mu_star)
return C_mu, regret
def Exploitation_MAB(MAB, vac_results, m=M, n=N):
arg_max = getBestResult(vac_results)
best_arm = getBestVaccin(MAB)
Ti = n // len(MAB)
regrets = []
rewards = (vac_results[arg_max] * Ti)
for m_i in range(m):
vaccins = MAB[arg_max].sample(1)
rewards += vaccins
regret = (n + m_i) * MAB[best_arm].immunity_rate - rewards
regrets.append(regret.item())
return regrets
def gloutonne_MAB(MAB, n=N, m=M):
exploration = Exploration(MAB, N)
regrets = Exploitation_MAB(MAB, exploration, m=m, n=n)
return regrets
#Comparaison
def compare():
MAB = generate_arms(K)
regrets_lay_robbins = []
C_mus = []
for i in range(1, M+1):
C_mu, regret_lay_robbins = lai_robbins(MAB, i)
C_mus.append(C_mu)
regrets_lay_robbins.append(regret_lay_robbins)
regrets_gloutonne = gloutonne_MAB(MAB, n=N, m=M)
return C_mus, regrets_lay_robbins, regrets_gloutonne
def plot_compare(C_mus, regret_lay_robbins, regrets_gloutonne):
plt.figure(figsize=(10, 6))
plt.plot((range(1, M+1)), C_mus, linestyle='-', label='C_mu')
plt.plot((range(1, M+1)), regret_lay_robbins, linestyle='-', label='Regret Lay & Robbins')
plt.plot((range(1, M+1)), regrets_gloutonne, linestyle='-', label='Regret Gloutonne')
plt.legend()
plt.title('Comparaison du regret Lay & Robbins et Gloutonne')
plt.grid(True)
plt.autoscale()
plt.savefig("figures/Lay_Robbins.png")
plt.show()C_mus, regrets_lay_robbins, regrets_gloutonne = compare()print(f'Moyenne du regret de l\'algorithme glouton: {np.mean(regrets_gloutonne)}')
print(f'Moyenne du regret de l\'algorithme Lay & Robbins: {np.mean(regrets_lay_robbins)}')
plot_compare(C_mus, regrets_lay_robbins, regrets_gloutonne)Moyenne du regret de l'algorithme glouton: 103.09383436203002
Moyenne du regret de l'algorithme Lay & Robbins: 91.20504772161748
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[13], line 3
1 print(f'Moyenne du regret de l\'algorithme glouton: {np.mean(regrets_gloutonne)}')
2 print(f'Moyenne du regret de l\'algorithme Lay & Robbins: {np.mean(regrets_lay_robbins)}')
----> 3 plot_compare(C_mus, regrets_lay_robbins, regrets_gloutonne)
Cell In[11], line 66, in plot_compare(C_mus, regret_lay_robbins, regrets_gloutonne)
63 plt.grid(True)
64 plt.autoscale()
---> 66 plt.savefig("figures/Lay_Robbins.png")
67 plt.show()
File /opt/conda/lib/python3.10/site-packages/matplotlib/pyplot.py:1023, in savefig(*args, **kwargs)
1020 @_copy_docstring_and_deprecators(Figure.savefig)
1021 def savefig(*args, **kwargs):
1022 fig = gcf()
-> 1023 res = fig.savefig(*args, **kwargs)
1024 fig.canvas.draw_idle() # Need this if 'transparent=True', to reset colors.
1025 return res
File /opt/conda/lib/python3.10/site-packages/matplotlib/figure.py:3378, in Figure.savefig(self, fname, transparent, **kwargs)
3374 for ax in self.axes:
3375 stack.enter_context(
3376 ax.patch._cm_set(facecolor='none', edgecolor='none'))
-> 3378 self.canvas.print_figure(fname, **kwargs)
File /opt/conda/lib/python3.10/site-packages/matplotlib/backend_bases.py:2366, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2362 try:
2363 # _get_renderer may change the figure dpi (as vector formats
2364 # force the figure dpi to 72), so we need to set it again here.
2365 with cbook._setattr_cm(self.figure, dpi=dpi):
-> 2366 result = print_method(
2367 filename,
2368 facecolor=facecolor,
2369 edgecolor=edgecolor,
2370 orientation=orientation,
2371 bbox_inches_restore=_bbox_inches_restore,
2372 **kwargs)
2373 finally:
2374 if bbox_inches and restore_bbox:
File /opt/conda/lib/python3.10/site-packages/matplotlib/backend_bases.py:2232, in FigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>(*args, **kwargs)
2228 optional_kws = { # Passed by print_figure for other renderers.
2229 "dpi", "facecolor", "edgecolor", "orientation",
2230 "bbox_inches_restore"}
2231 skip = optional_kws - {*inspect.signature(meth).parameters}
-> 2232 print_method = functools.wraps(meth)(lambda *args, **kwargs: meth(
2233 *args, **{k: v for k, v in kwargs.items() if k not in skip}))
2234 else: # Let third-parties do as they see fit.
2235 print_method = meth
File /opt/conda/lib/python3.10/site-packages/matplotlib/backends/backend_agg.py:509, in FigureCanvasAgg.print_png(self, filename_or_obj, metadata, pil_kwargs)
462 def print_png(self, filename_or_obj, *, metadata=None, pil_kwargs=None):
463 """
464 Write the figure to a PNG file.
465
(...)
507 *metadata*, including the default 'Software' key.
508 """
--> 509 self._print_pil(filename_or_obj, "png", pil_kwargs, metadata)
File /opt/conda/lib/python3.10/site-packages/matplotlib/backends/backend_agg.py:458, in FigureCanvasAgg._print_pil(self, filename_or_obj, fmt, pil_kwargs, metadata)
453 """
454 Draw the canvas, then save it using `.image.imsave` (to which
455 *pil_kwargs* and *metadata* are forwarded).
456 """
457 FigureCanvasAgg.draw(self)
--> 458 mpl.image.imsave(
459 filename_or_obj, self.buffer_rgba(), format=fmt, origin="upper",
460 dpi=self.figure.dpi, metadata=metadata, pil_kwargs=pil_kwargs)
File /opt/conda/lib/python3.10/site-packages/matplotlib/image.py:1689, in imsave(fname, arr, vmin, vmax, cmap, format, origin, dpi, metadata, pil_kwargs)
1687 pil_kwargs.setdefault("format", format)
1688 pil_kwargs.setdefault("dpi", (dpi, dpi))
-> 1689 image.save(fname, **pil_kwargs)
File /opt/conda/lib/python3.10/site-packages/PIL/Image.py:2429, in Image.save(self, fp, format, **params)
2427 fp = builtins.open(filename, "r+b")
2428 else:
-> 2429 fp = builtins.open(filename, "w+b")
2431 try:
2432 save_handler(self, fp, filename)
FileNotFoundError: [Errno 2] No such file or directory: 'figures/Lay_Robbins.png'

Réponse
Nous constatons que l'algorithme glouton obtient de meilleur résultat
que l'algorithme Lay & Robbins.
Cet algorithme améliore la version précédente en ajoutant un biais lié à la fréquentation de chaque vaccin :
$$ \bar{\mu}_i = \hat{\mu}_i + \sqrt{\frac{C\log{n}}{T_i}} $$,
avec C = 2.
Q8. Implémentez la modification de cet algorithme. Observez un intérêt à conserver N > 0 ? Et ϵ < 1 ? Expliquez pourquoi.
Dans la suite, on prendra N = 0 et ϵ = 1.
def ApplyBias_GetBestVaccin(vac_results, hist, n, C=2):
rewards_b = np.zeros(K)
for i in range(K):
if (hist[i] != 0):
rewards_b[i] = vac_results[i] / hist[i] + math.sqrt(C * math.log(n) / hist[i])
else:
rewards_b[i] = np.inf
return np.argmax(rewards_b)
def Exploitation_UCB(MAB, immunity_rate, n, m, eps, C=2):
Ti = n // K
hist = [Ti] * K
rewards = [immunity_rate[i] * hist[i] for i in range(K)]
for m_i in range(1, m+1):
if (torch.rand(1) > eps):
best_arm = torch.randint(0, len(MAB), (1,)).item()
else:
best_arm = ApplyBias_GetBestVaccin(rewards, hist, n+m_i, C=C)
arm = MAB[best_arm].sample(1).item()
hist[best_arm] += 1
rewards[best_arm] = rewards[best_arm] + arm
return np.sum(rewards)
def UCB(MABs, n, m, eps, C=2):
regret = []
for MAB in MABs:
#MAB = generate_arms(K)
if (n == 0):
exploration = [0] * len(MAB)
else:
exploration = Exploration(MAB, n)
best_arm = getBestVaccin(MAB)
exploitation = Exploitation_UCB(MAB, exploration, n, m, eps, C=C)
r_N = (n + m) * MAB[best_arm].immunity_rate - exploitation
regret.append(r_N)
return regretdef multi_plot_regret(total_regrets):
colors = ['b', 'r', 'g', 'y']
labels = ['epsilon = 0.2', 'epsilon = 0.5', 'epsilon = 0.8', 'epsilon = 1']
plt.figure(figsize=(15, 8))
for i in range(len(total_regrets)):
plt.plot(np.arange(0, 200, K),total_regrets[i], marker='x', linestyle='--', color=colors[i], label=labels[i])
plt.legend()
plt.xlabel('Taille de N')
plt.ylabel('Regret')
plt.title('Courbe du Regret en fonction de N et epsilon')
plt.grid(True)
plt.show()def UCB_regret_test():
eps = [0.2, 0.5, 0.8, 1]
regrets = [[] for i in range(len(eps))]
n_total = N + M
MABs = [generate_arms(K) for i in range(100)]
for n in range(0, 200, K):
for i in range(len(eps)):
regret = UCB(MABs, n, n_total-n, eps[i], C=2)
regrets[i].append(np.mean(regret))
return regrets
UCB_regret = UCB_regret_test()multi_plot_regret(UCB_regret)
Réponse
N > 0 indique la taille du training set N (nombre de vaccins dans la
phase d'exploration).
Ici, en appliquant l'algorithme UCB, les vacins les moins représentés
sont forcés à être testé, ce qui rend inutile la phase d'exploration
(les informations des différents vaccins sont constamment mis à
jour).
ϵ < 1 indique que lors de
la phase d'exploitation, il y a une chance ϵ de prendre un vaccin
aléatoirement, ce qui fait que les différents vaccins peuvent être tous
être sélectionnés.
Cependant, comme expliqué, l'algorithme UCB applique déjà ces procédés.
N > 0 et ϵ < 1 n'ont pas
d'intéret à être conservés.
Q9. Tracez sous la forme d'une animation l'évolution des taux d'immunisation empirique (fig. de gauche) et l'évolution du regret (fig. droite). Dans la figure de gauche, vous representerez μ̄i et μ̂i pour chaque vaccin.
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
import torch
def ApplyBias_GetBestVaccin_figure(vac_results, hist, n, C=2):
rewards_b = []
for i in range(len(vac_results)):
if hist[i] != 0:
bias = vac_results[i] + math.sqrt(C * math.log(n) / hist[i])
else:
bias = vac_results[i]
rewards_b.append(bias)
index = rewards_b.index(max(rewards_b))
return rewards_b, index
def Exploitation_proba_figure_animation(MAB, exploration, n, m, eps, C=2):
vaccins = []
total_immunity_rate = []
total_immunity_rate_bias = []
immunity_rate = exploration.copy() # Use a copy to avoid references
regret = []
Ti = n // len(MAB)
hist = [Ti] * len(MAB)
best_vaccin = getBestVaccin(MAB)
for m_i in range(m):
total_immunity_rate.append(immunity_rate.copy())
if torch.rand(1).item() > eps:
best_arm = torch.randint(0, len(MAB), (1,)).item()
total_immunity_rate_bias.append(immunity_rate.copy())
else:
rewards, best_arm = ApplyBias_GetBestVaccin_figure(immunity_rate, hist, n + m_i, C=C)
total_immunity_rate_bias.append(rewards)
arm = MAB[best_arm].sample(1)[0] # Removed .item(), get first element
hist[best_arm] += 1
immunity_rate[best_arm] = (immunity_rate[best_arm] * (hist[best_arm] - 1) + arm) / hist[best_arm]
vaccins.append(arm)
r_N = (n + m_i) * MAB[best_vaccin].immunity_rate - ((exploration[best_arm] * Ti) + sum(vaccins))
regret.append(r_N)
vaccins = [1 if vaccin else 0 for vaccin in vaccins]
return regret, total_immunity_rate, total_immunity_rate_bias
def UCB_figure_animation(MABs, n, m, eps, C=2):
all_regret = []
all_total_immunity_rate = []
all_total_immunity_rate_bias = []
for MAB in MABs:
if n == 0:
exploration = [0.0] * len(MAB)
else:
exploration = Exploration(MAB, n)
regret, total_immunity_rate, total_immunity_rate_bias = Exploitation_proba_figure_animation(MAB, exploration, n, m, eps, C=C)
all_regret.append(regret)
all_total_immunity_rate.append(total_immunity_rate)
all_total_immunity_rate_bias.append(total_immunity_rate_bias)
return all_regret, all_total_immunity_rate, all_total_immunity_rate_bias
# Animation Setup
def create_animation():
K = 5 # Number of vaccines
M = 100 # Total time
N = 50 # Initial time
eps = 0.1 # Exploration rate
C = 2 # Confidence parameter
# Generate arms
MABs = [generate_arms(K) for _ in range(1)]
# Calculate regrets and immunity rates
n_init = 5
m_init = M + N - n_init
regrets, total_immunity_rate, total_immunity_rate_bias = UCB_figure_animation(MABs, n_init, m_init, eps, C=C)
# Extract results (since MABs contains only one list)
regrets = regrets[0]
total_immunity_rate = total_immunity_rate[0]
total_immunity_rate_bias = total_immunity_rate_bias[0]
# Create subplots
fig, (ax_left, ax_right) = plt.subplots(1, 2, figsize=(16, 6))
# Left plot: Empirical immunization rates
x = np.arange(K)
width = 0.35 # Width of the bars
# Initialize bars
bars_mu = list(ax_left.bar(x - width/2, [0]*K, width, label=r'$\bar{\mu}_i$'))
bars_hat_mu = list(ax_left.bar(x + width/2, [0]*K, width, label=r'$\hat{\mu}_i$'))
ax_left.set_xlim(-0.5, K-0.5)
ax_left.set_ylim(0, 1) # Adjust based on immunity rates
ax_left.set_title("Évolution des taux d'immunisation empirique")
ax_left.set_xlabel("Vaccin")
ax_left.set_ylabel("Taux d'immunisation")
ax_left.set_xticks(x)
ax_left.set_xticklabels([f'V{k+1}' for k in range(K)])
ax_left.legend()
# Right plot: Regret curve
line_regret, = ax_right.plot([], [], color='red', label='Regret')
ax_right.set_title("Évolution du regret")
ax_right.set_xlabel("Temps")
ax_right.set_ylabel("Regret")
ax_right.set_xlim(0, len(regrets))
ax_right.set_ylim(0, max(regrets) * 1.1 if len(regrets) > 0 else 1) # Add margin, handle empty
ax_right.legend()
# Update function for animation
def update(frame):
if frame == 0:
return bars_mu + bars_hat_mu + [line_regret]
# Update bar heights
current_mu = total_immunity_rate[:frame]
mean_mu = np.mean(current_mu, axis=0)
current_mu_bias = total_immunity_rate_bias[:frame]
mean_mu_bias = np.mean(current_mu_bias, axis=0)
for k in range(K):
bars_mu[k].set_height(mean_mu[k])
bars_hat_mu[k].set_height(mean_mu_bias[k])
# Update regret line
x_regret = np.arange(frame)
y_regret = regrets[:frame]
line_regret.set_data(x_regret, y_regret)
return bars_mu + bars_hat_mu + [line_regret]
# Create animation
ani = FuncAnimation(fig, update, frames=len(regrets)+1, blit=True, interval=100, repeat=False)
# Display animation
return HTML(ani.to_jshtml())
# Run Animation
animation_html = create_animation()
display(animation_html)
Q10. Reprenez la question Q5 avec cette algorithme. Concluez sur l'utilité (ou l'inutilité) de la phase d'exploration. Comparez les performances d'UCB avec celles de l'algorithme glouton.
def UCB_regret_N_M():
eps = 1
regrets_total = []
MABs = [generate_arms(K) for i in range(1)]
n = K
m = N+M
while n <= m:
regret = UCB(MABs, n, m, eps, C=2)
n += K
m -= K
regrets_total.append(regret)
return regrets_total
regrets_UCB_SizeN = UCB_regret_N_M()
plot_regret(regrets_UCB_SizeN)
Réponse
N > 0 indique la taille du training set N (nombre de vaccins dans la
phase d'exploration).
Ici, en appliquant l'algorithme UCB, les vacins les moins représentés
sont forcés à être testé, ce qui rend inutile la phase d'exploration
(les informations des différents vaccins sont constamment mis à
jour).
La phase d"exploration est donc inutile
Q11. Testez différentes valeurs pour C et trouvez sa valeur optimale expérimentalement.
def UCB_C():
eps = 1
regrets_total = []
MABs = [generate_arms(K) for i in range(100)]
Cs = [0, 0.1, 0.5, 1, 2, 5, 10]
for c in Cs:
regret = UCB(MABs, 0, N+M, eps, C=c)
regrets_total.append(np.mean(regret))
min_index = regrets_total.index(min(regrets_total))
print (f'La meilleur valeur de C pour N = 0 et M={N+M}, est de C={Cs[min_index]} avec un regret moyen de {regrets_total[min_index]}')
UCB_C()La meilleur valeur de C pour N = 0 et M=550, est de C=0.1 avec un regret moyen de 14.630263328552246
Cet algorithme propose de modéliser la variable aléatoire de chaque vaccin avec une loi β dont les paramètres a et b correspondent au nombre de patients que le vaccin a immunisés (resp. non immunisés).
Pour chaque patient, on tire un valeur aléatoire pour la loi β décrivant chaque vaccin, puis on choisit le vaccin avec la plus grande valeur tirée.
Q12. Implémentez cet algorithme. En testant plusieurs valeurs de N, montrez que la phase d'exploration précédente a un impact très limité. Cela veut-il dire que l'algorithme ne contient pas d'initialisation ?
from torch.distributions import Beta
def Thompson_sampling(MAB, n):
# Initialisation
alpha = torch.ones(K)
beta = torch.ones(K)
rewards = []
# Phase d'exploration
for _ in range(n):
samples = torch.tensor([Beta(alpha[i], beta[i]).sample() for i in range(K)])
chosen_arm = torch.argmax(samples).item()
reward = MAB[chosen_arm].sample(1).item()
# Mise à jour des paramètres α et β pour le bras sélectionné
alpha[chosen_arm] += reward
beta[chosen_arm] += 1 - reward
rewards.append(reward)
# Phase d'exploitation
for _ in range(M):
samples = torch.tensor([Beta(alpha[i], beta[i]).sample() for i in range(K)])
chosen_arm = torch.argmax(samples).item()
reward = MAB[chosen_arm].sample(1).item()
# Mise à jour des paramètres α et β pour le bras sélectionné
alpha[chosen_arm] += reward
beta[chosen_arm] += 1 - reward
rewards.append(reward)
# Calcul du regret
best_arm = getBestVaccin(MAB)
regret = (n + M) * MAB[best_arm].immunity_rate - sum(rewards)
return regret
# Test de l'algorithme sur plusieurs instances
def test_Thompson_sampling(n_values=[0, 20, 50, 100]):
regrets = {}
for n in n_values:
regrets_n = []
for _ in range(number_MAB):
MAB = generate_arms(K)
regret = Thompson_sampling(MAB, n)
regrets_n.append(regret)
regrets[n] = regrets_n
return regrets
# Tester avec différentes valeurs de N
regrets_Thompson = test_Thompson_sampling(n_values=[0, 20, 50, 100])
# Comparaison des performances
for n, regret in regrets_Thompson.items():
print(f"N = {n}, Moyenne du regret: {np.mean(regret)}, Écart-type du regret: {np.std(regret)}")N = 0, Moyenne du regret: 14.319062232971191, Écart-type du regret: 11.961745262145996
N = 20, Moyenne du regret: 14.019129753112793, Écart-type du regret: 10.967537879943848
N = 50, Moyenne du regret: 13.356916427612305, Écart-type du regret: 11.429512023925781
N = 100, Moyenne du regret: 13.294836044311523, Écart-type du regret: 12.096404075622559
Q13. Tracez sous la forme d'une animation l'évolution des taux d'immunisation empirique (fig. de gauche) et l'évolution du regret (fig. droite). Dans la figure de gauche, vous representerez le taux d'immunisation empirique pour chaque vaccin avec un graphique en violon qui représente la loi beta associée à chaque vaccin.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from torch.distributions import Beta
K = 5
M = 500
N = 50
def Thompson_sampling_with_tracking(MAB, n):
alpha = torch.ones(K)
beta = torch.ones(K)
rewards = []
regret_history = []
mu_hat_history = []
best_arm = getBestVaccin(MAB)
for _ in range(n):
samples = torch.tensor([Beta(alpha[i], beta[i]).sample() for i in range(K)])
chosen_arm = torch.argmax(samples).item()
reward = MAB[chosen_arm].sample(1).item()
alpha[chosen_arm] += reward
beta[chosen_arm] += 1 - reward
rewards.append(reward)
mu_hat = alpha / (alpha + beta)
mu_hat_history.append(mu_hat.clone().detach())
regret = (n + M) * MAB[best_arm].immunity_rate - sum(rewards)
regret_history.append(regret)
for _ in range(M):
samples = torch.tensor([Beta(alpha[i], beta[i]).sample() for i in range(K)])
chosen_arm = torch.argmax(samples).item()
reward = MAB[chosen_arm].sample(1).item()
alpha[chosen_arm] += reward
beta[chosen_arm] += 1 - reward
rewards.append(reward)
mu_hat = alpha / (alpha + beta)
mu_hat_history.append(mu_hat.clone().detach())
regret = (n + M) * MAB[best_arm].immunity_rate - sum(rewards)
regret_history.append(regret)
return alpha, beta, mu_hat_history, regret_history
MAB = generate_arms(K)
alpha, beta, mu_hat_history, regret_history = Thompson_sampling_with_tracking(MAB, N)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
ax1.set_title("Taux d'immunisation empirique (\u03BC_i) et \u03B2(\u03B1_i, \u03B2_i)")
ax1.set_xlabel("Vaccins")
ax1.set_ylabel("Taux d'immunisation")
ax2.set_title("Évolution du regret")
ax2.set_xlabel("Itération")
ax2.set_ylabel("Regret")
def update(frame):
ax1.clear()
ax2.clear()
data = [Beta(alpha[i], beta[i]).sample((1000,)).numpy() for i in range(K)]
ax1.violinplot(data, showmeans=True, showmedians=True)
ax1.plot(np.arange(1, K + 1), mu_hat_history[frame], 'r--o', label="Taux empiriques (\u03BC_i)")
ax1.set_ylim(0, 1)
ax1.legend()
ax2.plot(regret_history[:frame], 'b-')
ax2.set_xlim(0, len(regret_history))
ax2.set_ylim(0, max(regret_history))
ani = FuncAnimation(fig, update, frames=len(mu_hat_history), interval=200, repeat=False)
# Peut prendre un peu de temps
display(HTML(ani.to_jshtml()))